# SETWD ===========================
rm(list=ls())


# LIBRARIES ===========================

library(tidyverse)
library(foreign)
library(stargazer)
library(lmtest)
library(multiwayvcov)
library(logistf)


# IMPORT DATA ===========================

dfgi <- read.csv("data/sve-08_04_24.csv",stringsAsFactors = F)


# ELECTION YEARS ANALYSIS ===========================

temp <- dfgi%>%
  filter(!is.na(gini_disp))%>%
  filter(elec.any==1)

## Descriptive stats ----
# percent of observations eroding:
mean(temp$erosion.strict)

## TAB A1. Logit: Election Years ----

m1 <- glm(erosion.strict ~ gini_disp, data=temp, family=binomial)
crse1 <- coeftest(m1, cluster.vcov(m1, temp$country.name))
crse1

m2 <- glm(erosion.strict ~ gini_disp + log(gdppc), data=temp, family=binomial)
crse2 <- coeftest(m2, cluster.vcov(m2, temp$country.name))
crse2

m3 <- glm(erosion.strict ~ gini_disp + log(gdppc) + year, 
          data=temp, family=binomial)
crse3 <- coeftest(m3, cluster.vcov(m3, temp$country.name))
crse3


stargazer(crse1,crse2,crse3, title="Logit: Election Years",digits=3,
          header=F,df=F,omit.stat=c("f","ser","aic","ll"),
          add.lines = list(c("Observations", nobs(crse1), nobs(crse2), nobs(crse3))),
          star.cutoffs=c(0.05,0.01,0.001),
          star.char=c("*","**","***"),
          notes="$^{*} p<0.05$; $^{**} p<0.01$; $^{***} p<0.001$",
          covariate.labels=c("Gini","Logged GDP per capita",
                             "Year"),
          notes.append=F,table.placement="h!",
          out="tables/apptab-elecs-10_27_24.tex")



# RARE EVENTS LOGIT ===========================

## TAB A2. Rare Events Logit ----

temp <- dfgi%>%
  filter(!is.na(gini_disp))

m1 <- logistf::logistf(erosion.strict ~ gini_disp, data = temp,
                       flic = TRUE)
summary(m1)

temp <- dfgi%>%
  filter(!is.na(gdppc),
         !is.na(gini_disp))

m2 <- logistf::logistf(erosion.strict ~ gini_disp + log(gdppc), data = temp,
                       flic = TRUE)
summary(m2)


m3 <- logistf::logistf(erosion.strict ~ gini_disp + log(gdppc) + year, data = temp,
                       flic = TRUE, control=logistf.control(maxit=1000))

summary(m3)


## TAB A3. Rare Events Logit: Election Years ----

temp <- dfgi%>%
  filter(!is.na(gini_disp)&elec.any==1)

m1 <- logistf::logistf(erosion.strict ~ gini_disp, data = temp,
                       flic = TRUE)
summary(m1)

temp <- dfgi%>%
  filter(!is.na(gdppc)&
           !is.na(gini_disp)&elec.any==1)

m2 <- logistf::logistf(erosion.strict ~ gini_disp + log(gdppc), data = temp,
                       flic = TRUE)
summary(m2)


m3 <- logistf::logistf(erosion.strict ~ gini_disp + log(gdppc) + year, data = temp,
                       flic = TRUE, control=logistf.control(maxit=1000))

summary(m3)



## TAB A4. Haggard and Kaufman Recode ----

temp <- dfgi%>%
  filter(!is.na(gini_disp))

temp2 <- temp
temp2$erosion.strict[temp2$country.name%in%c("Philippines","Senegal","India")] <- 0

m.exc <- glm(erosion.strict ~ gini_disp, data=temp2, family=binomial)
mhk1 <- coeftest(m.exc, cluster.vcov(m.exc, temp2$country.name))
mhk1

m.exc <- glm(erosion.strict ~ gini_disp + log(gdppc), data=temp2, family=binomial)
mhk2 <- coeftest(m.exc, cluster.vcov(m.exc, temp2$country.name))
mhk2

m.exc <- glm(erosion.strict ~ gini_disp + log(gdppc) + year, data=temp2, family=binomial)
mhk3 <- coeftest(m.exc, cluster.vcov(m.exc, temp2$country.name))
mhk3

stargazer(mhk1,mhk2,mhk3, title="Haggard and Kaufman Recode",digits=3,
          header=F,df=F,omit.stat=c("f","ser","aic","ll"),
          add.lines = list(c("Observations", nobs(mhk1), nobs(mhk2), nobs(mhk3))),
          star.cutoffs=c(0.05,0.01,0.001),
          star.char=c("*","**","***"),
          notes="$^{*} p<0.05$; $^{**} p<0.01$; $^{***} p<0.001$",
          covariate.labels=c("Gini","Log(GDPpc)","Year"),
          notes.append=F,table.placement="h!",
          out="tables/apptab-hk-10_27_24.tex")


